A damage-plasticity model for the dynamic failure of concrete 
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ABSTRACT: A constitutive model based on the combination of damage mechanics and plasticity is developed to analyse concrete 
structures subjected to dynamic loading. The aim is to obtain a model, which requires input parameters with clear physical 
meanings. The model should describe the important characteristics of concrete subjected to multiaxial and rate-depending loading. 
This is achieved by combining an effective stress based plasticity model with an isotropic damage model based on plastic and 
elastic strain measures. The model response in tension, uni-, bi- and tri-axial compression is compared to experimental results in 
the literature. 
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1 INTRODUCTION 

Concrete is a strongly heterogeneous material, which exhibits a 
complex nonlinear mechanical behaviour Failure in tension and 
low confined compression is characterised by softening which 
is defined as decreasing stress with increasing deformations. 
This softening response is accompanied by a reduction of the 
unloading stiffness of concrete, and irreversible (permanent) 
deformations, which are localised in narrow zones often 
called cracks or shear bands. On the other hand, the 
behaviour of concrete subjected to high confined compression is 
characterised by a ductile hardening response; that is, increasing 
stress with increasing deformations. Furthermore, high loading 
rates are known to significantly increase the strength in tension 
and compression. These phenomena should be considered in 
a constitutive model for analysing the dynamic behaviour of 
concrete structures. 

There are many constitutive models for the nonlinear 
response of concrete proposed in the literature. Commonly used 
frameworks are plasticity, damage mechanics and combinations 
of plasticity and damage mechanics. Stress-based plasticity 
models are useful for the modelling of concrete subjected to 
triaxial stress states, since the yield surface corresponds at a 
certain stage of hardening to the strength envelope of concrete. 
Furthermore, the strain split into elastic and plastic parts 
represents realistically the observed deformations in confined 
compression, so that unloading and path-dependency can be 
described well. However, plasticity models are not able to 
describe the reduction of the unloading stiffness that is observed 
in experiments. Conversely, strain based damage mechanics 
models are based on the concept of a gradual reduction of the 
elastic stiffness driven by strain measures. For isotropic damage 
mechanics models, the stress evaluation procedure is explicit, 
which allows for a direct determination of the stress state, 
without an iterative calculation procedure. Furthermore, the 



stiffness degradation in tensile and low confined compressive 
loading observed in experiments can be described. However, 
isotropic damage mechanics models are often unable to describe 
irreversible deformations observed in experiments and are 
mainly limited to tensile and low confined compression stress 
states ||T|. On the other hand, combinations of isotropic damage 
and plasticity are widely used and many different models have 
been proposed in the literature. One popular class of models 
relies on a combination of stress-based plasticity formulated in 
the effective stress space combined with a strain based damage 
model 12J, |i3J, \5j, ^7J. 



In the present work, the combined damage-plasticity model 
presented in (5), (8) (CDPMl) is revisited to develop an 
constitutive model for the rate dependent failure of concrete, 
which is characterised by its numerical stability, well defined 
input parameters and flexibility to be adapted to newly 
developed concrete based materials, such as fibre reinforced 
concrete. As this model can be seen as an augmentation of 
CDPMl it is called here CDPM2. The stress-based plasticity 
part of the model is based on the effective stress. It includes 
hardening in the post-peak regime, which is used to model 
the strain rate dependence of strength by delaying the onset 
of damage. The plasticity part is combined with a damage 
model, which is based on elastic and plastic strain measures 
and distinguishes between tensile and compressive stress states 
using an approach similar to those proposed in yj, [^J, [lOj . The 
damage model is used to describe the complex strength envelope 
of concrete. With this combination of plasticity and damage 
mechanics, it is aimed to provide a computationally efficient 
model for the dynamic behaviour of concrete. 



2 MODEL 

2.1 General framework 

The damage plasticity constitutive model is based on the 
following stress-strain relationship: 

f7 = (1 -C(Jt)at + (l -«c)^c (1) 

where Oi and are the positive and negative parts of the 
effective stress tensor a, respectively, and (Oi and a)c are two 
scalar damage parameters, ranging form (undamaged) to 1 
(fully damaged). The effective stress a is defined as 



a, a 



e-en 



(2) 



where Dg is the elastic stiffness tensor based on the elastic 
Young's modulus E and Poisson's ratio v, £ is the strain tensor 
and Ep is the plastic strain tensor The positive and negative 
parts of the effective stress a in Eq. ([T]) are determined from 
the principal effective stress as ffpt = (o'p)^ and dpc = (ffp) , 
where ()^ and ()_ are positive and negative part operators, 
respectively, defined as (jc)^ — max (0,x) and = min (0,ji:). 
In addition, a scalar measure is defined which distinguishes 
between tensile and compressive stress states. It has the form 



(ap,-)_((ap,)^ + (ap,)_) 



(3) 



where (o'pi)^^ and (o'pi) are the components of the compres- 
sive and tensile part of the principal stresses, respectively. The 
parameter Uq varies in the range from to 1 . For instance, for 
a combined tensile and compressive stress state with principal 
stress components ffp = (— d',0.2d',0.1d')^, the positive and 
negative principal stresses are apt = (0,0.2d',0.1d')^ and ffpc — 
(—a, 0,0)^, respectively. For this stress state, the variable is 
Uc = 0.95. This measure is used later in the definition of the 
damage parameter. 

The plasticity model is based on the effective stress and thus 
independent of damage. The model is described by the yield 
function, the flow rule, the evolution law for the hardening 
variable and the loading-unloading conditions. The form of the 
yield function is 



/p ((T,K'p) =f ((J,^h,?s) 



(4) 



where and are hardening functions, which depend on the 
plastic hardening parameter K"p. The flow rule is 



^ da 



(5) 



where Cp is the rate of the plastic strain, X is the rate of the 
plastic multiplier and gp is the plastic potential. The rate of the 
hardening parameter fCp is related to the rate of the plastic strain 
by an evolution law, which is presented in Section 2.2 The 
loading-unloading conditions are 



/p < 0, A > 0, A/p = 



(6) 



A detailed description of the individual components of the 




Figure 1 . Geometrical meaning of the inelastic strain for the 
combined damage-plasticity model. The inelastic strain is 
composed of reversible O) (e — Cp) and irreversible Cp parts. 
The dashed lines represent elastic unloading with the same 
stiffness as the initial elastic loading. 



The damage part of the present damage-plasticity model is 
related to elastic and plastic strain measures. For a pure tensile 
stress state, fft = O, Oc is zero and a)t = ft), so that Eq. ([TJ 
becomes 



C7 = (1 - ftj) C7 = (1 - ft))De : (e -ep 



(7) 



The equation is rearranged as 



(J=De 



e- En 



ftj 



(e-ep)))=De:(e-ei) (8) 



where 

ei = Ep + ft) (e - Ep) (9) 

is the inelastic strain which is subtracted from the elastic strain. 
The geometrical interpretation of the inelastic strain and its 
split for uniaxial tension, hardening plasticity and linear damage 
evolution are shown in Fig.[T| The part ft) (e — Ep) is reversible 
and Ep is irreversible. These two parts of the inelastic strain are 
used to define the damage history variables, see Section 2.3 



plasticity model are discussed in Section 2.2 



2.2 Plasticity part 

The plasticity model is formulated in a three-dimensional 
framework with a pressure-sensitive yield surface, hardening 
and non-associated flow. The main components are the yield 
condition, the hardening law, the evolution law for the hardening 
variable and the flow rule. This model is an extension of 
CDPMl in l|5l, in which the plasticity response is assumed to 
be perfect-plastic in the regime in which damage is active. In 
the present study, the plasticity part exhibits hardening in this 
regime, which is used to model the strain rate dependence of 
the strength by adjusting the onset of damage. This extension to 
the hardening requires several extensions of the plasticity part 
of CDPMl, which is presented in the following section. 

The yield surface is described in terms of the cylindrical 
coordinates in the principal effective stress space (Haigh- 
Westergaard coordinates), which are the volumetric effective 




Figure 2. The evolution of the meridional section of the yield 
surface during hardening in the pre-peak and post-peak 
regime. 




Figure 3. The evolution of the deviatoric section of the yield 
surface during hardening for a constant volumetric stress 
of ay = -/c/3. 



stress ffy, the norm of the deviatoric stress p and the Lode angle 
of the deviatoric stress 9. For a definition of these coordinates it 
is referred to |5|. 
The yield function 



/p(ffv,p,0;K-p 



(10) 



^r(cos0) + ^ 

/6/c /e 



-?h(K"p)?s(K-p) 



depends on the effective stress (which enters in the form of 
cylindrical coordinates) and on the hardening variable JCp (which 
enters through the dimensionless variables and q^,). Parameter 
/c is the uniaxial compressive strength. 

The meridians of the yield surface /p = are parabolic, and 
the deviatoric sections change from triangular shapes at low 
confinement to almost circular shapes at high confinement. The 
shape of the deviatoric section is controlled by the function 



r(cos 9) = 



4(l-e2)cos2 + (2e-l)2 



(11) 



2(1 - e2) cos e + (2e - 1 ) V'4(l -e2)cos2 +5e2 -4e 

proposed by pT[ . The eccentricity parameter e and the friction 
parameter 

mo^ , ^ (12) 

/c/t e + 1 

are calibrated from the values of uniaxial and equibiaxial 
compressive strengths and uniaxial tensile strength as described 
in |12|. The shape of the meridians of the yield surface is 
controlled by the hardening variables ^j, and and mo. The 
evolution of the yield surface during hardening is presented 
in Figure |2] and |3] The two hardening functions, which are 
functions of the hardening variable fCp are presented in Figure |4] 
If the two variables and are equal to one, the yield surface 
turns into the failure surface proposed by [13 1. 

In the present model, the flow rule in Eq. (jSj) is non- 
associated, which means that the yield function /p and the 
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Figure 4. The two hardening laws (solid line) and q^ (dashed 
line). 



plastic potential g-p do not coincide and, therefore, the direction 
of the plastic flow m = dgp/da is not normal to the yield 
surface. This is important for realistic modeling of the 
volumetric expansion under compression for frictional materials 
such as concrete. An associated flow rule for this type of yield 
surface gives an unrealistically high volumetric expansion in 
compression, which leads in the case of passive confinement to 
an overestimated strength (peak stress); see p4) . The plastic 
potential is given as 



|[l-?h(K-p)] 
+ qh{Kp)qsiKp) 



%/6/c /c 



\V6fc 





/If 






mg((7v, 




/c ; 



(13) 



where 



'tg(CTy, K-p) = Ag (fC-p) Bg (fC-p) /eexp , ( 



(14) 



is a variable controlling the ratio of volumetric and deviatoric 
plastic flow. The variables Ag (fCp) and Bg (ffp), which depend 
on the hardening function q^, are derived from assumptions 
on the plastic flow in uniaxial tension and compression. In 
uniaxial tension, the laterial plastic strain is chosen to be zero. In 



uniaxial compression, the ratio of axial and laterial plastic strain 
is set to the model constant D{. The plastic potential does not 
depend on the third Haigh-Westergaard coordinate (Lode angle 
0). This increases the efficiency of the implementation and the 
robustness of the model. However, it also limits the capability of 
this flow rule to describe the response of concrete in multiaxial 
compression. 

The evolution law for the hardening variable Kp sets the rate of 
the hardening variable equal to the norm of the plastic strain rate 
scaled by a hardening ductility measure. This scaling factor is 
constructed such that the model response is more ductile under 
compression |5J. 

2.3 Damage part 

The damage parameter is determined by means of history 
variables that are based on measures of the plastic and elastic 
strain. The measure of the plastic strain is based on the norm 
of the rate of the plastic strain. The measure of the elastic 
strain, denoted here as the equivalent strain £, is more difficult 
to choose. With the equivalent strain, the onset of damage 
is determined. Thus, it is required to choose an equivalent 
strain measure, which represents at a certain value stress states 
which are located on the strength envelope of concrete. In 
the present work, the equivalent strain e is chosen so that the 
effective stress states at the onset of damage satisfy the strength 
criterion developed by Menetrey and Willam 1 13 1 on which the 
yield surface in the plasticity part is based. The corresponding 
damage envelope has the form 



above for uniaxial tension, the damage hardening parameter is 
determined as 

qd = ^ (19) 

/c 

Consequently, q^i is proportional to the effective stresses in 
tension and compression, respectively, and is equal to one, if 
the effective stresses are equal to the strength described by the 



criterion in Eq. ( 15 1. Therefore, the variable is well suited to 
be used for the definition of the equivalent strain, which is 



/t 



(20) 



Based on the plastic hardening parameter (Cp, and the 
equivalent strain e, six damage history variables (K^t, JQic, K'dit, 
"die. K'd2ti Kd2c) are defined. The evolution of these damage 
variables is controlled by two damage loading functions 



1 - 

Jdt = — £ - K^dt 
(Xi- 



and 



K'dc 



for which the loading-unloading conditions 



and 



/dt < 0, fc-dt > 0, K-dt/dt = 



fdc < 0, K-dc > 0, K-dc/dc = 



(21) 
(22) 

(23) 
(24) 



3p! 
2fi 



/6/c 



(COS0) 



CTv 

/c 



apply, respectively. With the variable in Eqs. (21 1 and (22 1 the 



(15) strain rate dependence of the material is modelled, see Sec. 2.4 



From Eq. (15i, the variable (positive root) is determined. 



which is used to define the equivalent strain e for the strain 
driven damage model. The meaning of q^ is illustrated by two 
representative effective stress states, namely uniaxial tension 
and uniaxial compression. For uniaxial tension, the effective 
stress state is defined as ffi = Oi, d'2 = O3 = 0, Oy ~ O't/3, 
si = 2at/3, S2 = S2 = — O't/3 and p = ^2/3at and cos = 1. 



Furthermore, from Eq. ([3]) is used to distinguish between 
tensile and compressive stress states. The next two variables 
are linked to the plastic hardening parameter. Here, only plastic 
strains after the onset of damage are considered. They are 
defined as 



Kdlt 



— K-p if KTdt > A K-dt > £0 

o' 



if K-dt = V K-dt < £0 



Setting this in Eq. ( 15 1 results in 



and 



/d(^t) 



' qs.m 




(16) 



1 



With r = - for = and the definition of otq in Eq. ( 12 1, this 

e 

simplifies to 



-K-p if K-dc > A K-dc > £0 
Kdlc = < J^s«r (26) 

[0 if K-dc = V K-dc < £0 

where jCg is a ductility measure, which describes the influence of 
multiaxial stress states on the softening response and controls so 
the ratio of fracture energies in tension Gft and compression G^. 
These two history variables represent the irreversible component 



/d(fft) 



fi 



- f - f 

-qd^i '2^' ~gd = o 

Jc/t 



The positive root of this equation for g'd is 

qd - 



/t 



(17) 



(18) 



of the inelastic strain introduced in Sec. 2. 1 Finally, the last two 
damage history variables are related to the maximum equivalent 
strains K-dt and K-dc and are defined as 



and 



On the other hand, for uniaxial compression, the effective 
stress state is defined as Gi = — Ob, 02 = 03= 0, ay = —0^/3 



K-d2t 



K-d2c 



Kdx 
K-dc 



(27) 



(28) 



and p = \/-(Jc- Repeating the steps, which were presented 



All six damage history variables are used to determine the two 
damage parameters Ot and cOc in Eq. ([T]i. The form of the two 




Wfi Wf 

Figure 5. Bilinear softening. 



damage parameters o, and Oc depends on the type of softening 
relation that is modelled. For instance, for linear softening in 
uniaxial tension, the stress inelastic displacement relationship 
in the softening regime is 







if < hSi < Wf 
if Wf < hei 



(29) 



where Wf is the crack opening at which the uniaxial stress is 
equal to zero. To extend this relationship to general loading, the 
inelastic strain Ci is replaced by the damage history variables as 



ei = K'dlt + ftJtK"d2t 



(30) 



Furthermore, the general stress-strain relationship in Eq. ([TJ 
reduces for uniaxial tension to 



a = {l-cOt)EKAt 



(31) 



where £ is replaced by JCdt to make it applicable to general strain 



states. Setting Eqs. (29i and (31 1 equal gives 



/tK-dit/i + Wf^K-dt - Wf/t 
K-dtWf£'-/tK-d2t/! 



(32) 



The expression for (Oc is obtained, if Kdt, K^n and Kd2t are 
replaced by fCdc, fQnc and K"d2c, respectively. Thus, 



CO, 



_ ftKdlch + WfE Kdc - Wf/t 
KdcWf£'-/tK-d2c/l 



(33) 



This procedure can be extended to derive the damage parameter 
for any combination of line segments, such as the bilinear 
stress crack opening curve shown in Figure [5] which is used 
in the analyses presented in this paper in Sections 2.5 and 3. 
For these analyses the parameters of the bilinear stress-strain 
curve were chosen as Ci — 0.2/t and Wfi — 0.2wf, for which 
Wf — 4.444GFt//t. For all the experiments, the axial specimen 
length was assumed to be ^^f = 0. 1 m for presenting the stress- 
strain curves. 

2.4 Rate dependence 

The response of concrete is strongly rate dependent. If the 
loading rate is increased, the apparent tensile and compressive 



strength increase. This increase is more pronounced in tension 
than in compression. In the present model, this rate dependence 



is taken into account in Eq. ( 23 i and ( 24 1 by the factor (X^ > 
1. The greater the rate factor ctr, the greater is the delay 
of the onset of damage and, therewith, the strength. For 
tension, the rate dependence is modeled using the expressions 
proposed by Malvar and Ross in |15|. For compression, the 
expressions reported in CEB-FIP Model Code 1990 |16| are 
used. Accordingly, the factor is defined as 



Or (1 -ac)art + acarc 



(34) 



where ttc is the factor introduced in Eq. (|3]l to model the different 
strain rate dependence in tension and compression. 
For tension, the rate factor is 



1 



art = < 



with 



and 



5s 



'max 

eto 



for gmax < 30 X 10-*^ s- 



for 30 X 10-*^ s-1 < e^ax < 1 s"' 



1/3 



for 1 s ' < a 



1 

l+8/c//c 



cO 



logjSs =65s-2 



(35) 
(36) 

(37) 



where taiax is the maximum principal strain rate component, 
/cO = 10 MPa and gto = 1 x 10"^ s^'. For uniaxial tension, 
for instance, the maximum principal strain rate component tmax 
is equal to the strain rate in uniaxial tension. For compression, 
the rate factor is 



1 



eco 



1.026as 



1/3 



for |emin, _ 
for 30 X 10-*^ s-i 

for30s-i < |e„i. 



< 30x 10-^ s-i 



< e,r 



where the parameter is defined as 

1 



0£s 



5 + 9/c//co 



and 



log 7s = 6.156as -2 



<30s 



(38) 



(39) 



(40) 



-1 



Here, is the minimum principal strain rate component and 
Eco = 30 X 10^^ s^'. For uniaxial compression, £n^n is equal to 
the axial compressive strain rate. 

2.5 Model response for varying strain rates and cyclic loading 

The response of the constitutive model is illustrated by several 
idealised load cases, before it is compared to a wide range 
of experimental results in the next section. Firstly, a quasi- 
static strain cycle is considered, which results in a stress-strain 
response shown in Fig. [6] The strain is increased from "0" to 
"1", where the tensile strength of the material is reached. Up 
to point "1", the material response is elastic-plastic with small 




plastic strains. With a further increase of the strain form "1" to 
"2", the effective stress part continues to increase, since //p > 0, 
whereas the nominal stress decreases, since the tensile damage 
variable COt increases. A reverse of the strain at point "2" results 
in an reduction of the stress with an unloading stiffness, which 
is less than the elastic stiffness of an elasto-plastic model, but 
greater than the stiffness of an elasto-damage mechanics model. 
At point "3"when the stress is equal to zero, a further reduction 
of the strain leads to a compressive response following a linear 
stress-strain relationship between the points "3" and "4" with the 
original Young's modulus E of the undamaged material. This 
change of stiffness is obtained by using two damage parameters, 
COt and (Bj.- At point "3" Ot > 0, but (0^ = 0. Up to "4" no 
farther plastic strains are generated, since the hardening from 
"0" to "1" has increased the elastic domain of the plasticity part 
so much that the yield surface is not reached at "4". A further 
decrease of the strain activates the compressive damage variable 
which leads to a reduction of the nominal stress. At point "5", 
the plasticity surface is reached and a subsequent increase of 
strain results in hardening of the plasticity part. However, the 
nominal stress, shown in Figure|6j decreases, since Oc increases. 
A second reversal of the strain direction ("6") changes the stress 
from tension to compression at "7", which is again associated 
with a change of the stiffness. 

The second group of examples consists of several tensile and 
compressive loading cases with constant strain strain rates. For 
uniaxial tension, strain rates of 1 x 10^^, 1, 10 and 100 1/s are 
applied. The corresponding stress strain responses are shown 
in Fig. |7] An increase of the loading rate results in a delay of 
the onset of damage. The strength is increased by factor (X^ in 
Eq. ( [34| , whereas the fracture energy is increased by the factor 
a^. The initial stiffness is modelled to be independent of the 
strain rate. 

For uniaxial compression, strain rates of — 1 x 10^^, —10 and 
— 100 are considered. The results are shown in Figure [8] Again, 
the strength in compression is increased by the factor and the 
fracture energy in compression by a^. In compression, the strain 
at peak stress for greater strain rates than for quasistatic loading 
depends strongly on the modulus Hp of the hardening function 

depicted in Fig. |4] An increase of the strain rate results in a 
delay of the onset of damage. This delay shifts the point of peak 
stress from intersection point of the hardening functions and 

into the hardening regime of q^^. 
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Figure 8. Rate effect in compression: Stress strain response for 
three compressive strain rates. 



3 COMPARISON WITH EXPERIMENTAL RESULTS 

In this section the model response is compared to five groups 
of experiments reported in the literature. For each group of 
experiments, the model constants E, v, /c, /t, Gpt and Gpc 
are adjusted to obtain a fit for the different types of concrete 
used in the experiments. The other model parameters are set 
to their default values. The hardening modulus (shown in 
Figure|4jl is set to = 0.5. The first analysis is uniaxial tensile 
setup with unloading. The model response is compared to the 
experimental results reported by Gopalaratnam and Shah in 1 17 1 
(Figure |9|. The relevant model constants for this experiment 
are £ = 25 GPa, v = 0.2, /c = 40 MPa, ft = 3.5 MPa, Gr = 
55 J/m2. 

The next example is an uniaxial compression test with 
unloading, for which the model response is compared to 
experimental results reported by Karsan and Jirsan [ISJ 
(Figure 10 1. The model constants are £ = 30 GPa, v = 0.2, 
/c = 28 MPa, /t = 2.8 MPa, Gpc = 2205 J/m^ . 

Next, the model is compared to uniaxial and biaxial 
compression tests reported by Kupfer et al. in |19|. For these 
experiments, the model constants are set to £ = 32 GPa, v = 0.2, 
/c = 32.8 MPa, /t = 3.3 MPa, Gr = 50 J/m^, Gpc = 3500 J/m^. 
The comparison with experimental results is shown in Figure 1 1 
for uniaxial, equibiaxial and biaxial compression. For the 
biaxial compression case, the stress ratio of the two compressive 
stress components is ai/G2 = —1/ — 0.5. 
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Figure 10. Uniaxial compression: Model response compared to Fig"^ ^3. hydrostatic compression: Model response compai-ed 
experimental results reported in [W] . to experiments used in |20) . 
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Figure 1 1 . Uniaxial compression: Model response compared to 
experimental results reported in p9) . 



Furthermore, the performance of the model is evaluated for 
triaxial tests and a hydrostatic test reported in [20 1. The 
material constants for this test are £ = 25 GPa, v = 0.2, 
/c = 28 MPa, /t = 2.8 MPa, Gpc = 2205 J/m^. The model 
response is compared to experimental results presented in 



Figures 12 and 13 



Finally, the model response in triaxial compression is 
compared to the experimental results reported in pT|. The 



material constants for this test are E — 25 GPa, v = 0.2, /c 
28 MPa, ft = 2.8 MPa, Gr 100 J/m^, Gpc = 15000 J/m^. 
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Figure 14. Confined compression: Model response compared 
to experiments reported in pTj . 



Overall, the agreement of the model response with the 
experimental results is very good. The model is able to represent 
the strength of concrete in tension and multiaxial compression. 
In addition, the strains at maximum stress in tension and 
compression agree well with the experimental results. The 
bilinear stress-crack opening curve that was used results in a 



good approximation of the softening curve in uniaxial tension 
and compression. 

4 CONCLUSIONS 

The present damage plasticity model combines a stress-based 
plasticity part with a strain based damage mechanics model. The 
model response is in good agreement with experimental results 
for a wide range of loading from uniaxial tension to confined 
compression. In the next steps, the model will be appUed to 
boundary value problems to evaluate its performance to describe 
failure processes of concrete mesh-independently. It is expected 
that the model will perform well for these problems since the 
softening response is formulated so that the crack band approach 
can be applied, which is known to lead in many cases to mesh- 
independent results. 
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